home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
EnigmA Amiga Run 1996 April
/
EnigmA AMIGA RUN 06 (1996)(G.R. Edizioni)(IT)[!][issue 1996-04][Skylink CD V].iso
/
earcd
/
utilgfx
/
raylab.lha
/
RayLab
/
source
/
intersct.c
< prev
next >
Wrap
C/C++ Source or Header
|
1996-01-04
|
13KB
|
422 lines
/*
name: intersct.c
Intersection-routines
---------------------
These routines will calculate the t-parameter of a line,
which corresponds to the point where the given line intersects with
the given primitive. The routines return zero (0) if no intersection
occured.
*/
#include <stdio.h>
#include "defs.h"
#include "extern.h"
/**********************************************************************
*
* Get (possible) intersection with given object
*
**********************************************************************/
double Intersect_Object(LINE *Line, OBJECT *Object)
{
LINE TransformedLine;
double t;
long i;
VECTOR tempv;
CopyPoint(&TransformedLine.Origin,&Line->Origin);
CopyVector(&TransformedLine.Direction,&Line->Direction);
if(Object->Transform.NumTransforms>0) {
for(i=Object->Transform.NumTransforms-1L;i>=0L;i--) {
switch(Object->Transform.Entry[i].Type) {
case TRANSFORM_SCALE:
TransformedLine.Origin.x=TransformedLine.Origin.x/Object->Transform.Entry[i].Values.x;
TransformedLine.Origin.y=TransformedLine.Origin.y/Object->Transform.Entry[i].Values.y;
TransformedLine.Origin.z=TransformedLine.Origin.z/Object->Transform.Entry[i].Values.z;
TransformedLine.Direction.x=TransformedLine.Direction.x/Object->Transform.Entry[i].Values.x;
TransformedLine.Direction.y=TransformedLine.Direction.y/Object->Transform.Entry[i].Values.y;
TransformedLine.Direction.z=TransformedLine.Direction.z/Object->Transform.Entry[i].Values.z;
break;
case TRANSFORM_MOVE:
TransformedLine.Origin.x=TransformedLine.Origin.x-Object->Transform.Entry[i].Values.x;
TransformedLine.Origin.y=TransformedLine.Origin.y-Object->Transform.Entry[i].Values.y;
TransformedLine.Origin.z=TransformedLine.Origin.z-Object->Transform.Entry[i].Values.z;
break;
case TRANSFORM_ROTATE:
NegVector(&tempv,&Object->Transform.Entry[i].Values);
RotatePoint(&TransformedLine.Origin,&tempv);
RotateVector(&TransformedLine.Direction,&tempv);
break;
case TRANSFORM_NONE:
default:
break;
}
}
}
switch(Object->ShapeType) {
case SHAPE_PLANE:
t=Intersect_Plane(&TransformedLine, (PLANE *)(Object->Shape));
break;
case SHAPE_SPHERE:
t=Intersect_Sphere(&TransformedLine, (SPHERE *)(Object->Shape));
break;
case SHAPE_ELLIPSOID:
t=Intersect_Ellipsoid(&TransformedLine, (ELLIPSOID *)(Object->Shape));
break;
case SHAPE_TRIANGLE:
t=Intersect_Triangle(&TransformedLine, (TRIANGLE *)(Object->Shape));
break;
case SHAPE_BOX:
t=Intersect_Box(&TransformedLine, (BOX *)(Object->Shape));
break;
case SHAPE_DISC:
t=Intersect_Disc(&TransformedLine, (DISC *)(Object->Shape));
break;
case SHAPE_CYLINDER:
t=Intersect_Cylinder(&TransformedLine, (CYLINDER *)(Object->Shape));
break;
default:
t=0.0;
}
return(t);
}
/**********************************************************************
*
* These are the actual intersection calculation routines
*
**********************************************************************/
double Intersect_Plane(LINE *Line, PLANE *Plane)
{
register double t,pnx,pny,pnz,a,tmp;
double lvx,lvy,lvz,lx0,ly0,lz0;
PlaneIntersectionAttempts++;
lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
pnx=Plane->Normal.x; pny=Plane->Normal.y; pnz=Plane->Normal.z; a=Plane->a;
tmp=pnx*lvx+pny*lvy+pnz*lvz;
if(tmp!=0.0) {
t=-(pnx*lx0+pny*ly0+pnz*lz0-a)/tmp;
if(t<EPSILON) { t=0.0; }
else { PlaneIntersections++; }
}
else { t=0.0; }
return(t);
}
double Intersect_Sphere(LINE *Line, SPHERE *Sphere)
{
double vx,vy,vz,lx0,ly0,lz0,sx0,sy0,sz0,r;
register double s,sroot,vxdx,vydy,vzdz,vxsqr,vysqr,vzsqr;
double t;
double dx,dy,dz;
double dxsqr,dysqr,dzsqr;
double vsum,t1,t2;
SphereIntersectionAttempts++;
vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
sx0=Sphere->Centre.x; sy0=Sphere->Centre.y; sz0=Sphere->Centre.z;
r=Sphere->r;
/* These are used more than once, */
dx=lx0-sx0; dy=ly0-sy0; dz=lz0-sz0; /* thus much time is saved by not */
vxdx=vx*dx; vydy=vy*dy; vzdz=vz*dz; /* performing the same operation */
dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz; /* twice (or even more). */
vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;
vsum=(vxsqr+vysqr+vzsqr);
s=-(vxdx+vydy+vzdz);
sroot=2*(vxdx*(vydy+vzdz)+vydy*vzdz)+r*r*vsum;
sroot-=vxsqr*(dysqr+dzsqr)+vysqr*(dxsqr+dzsqr)+vzsqr*(dxsqr+dysqr);
if(sroot>=0) {
sroot=sqrt(sroot);
t1=(s+sroot)/vsum;
t2=(s-sroot)/vsum;
if((t1>EPSILON)&&((t1<t2)||(t2<EPSILON))) {
t=t1;
SphereIntersections++;
}
else {
if((t2>EPSILON)&&((t2<t1)||(t1<EPSILON))) {
t=t2;
SphereIntersections++;
}
else { t=(double) 0; }
}
}
else t=(double) 0;
return(t);
}
double Intersect_Ellipsoid(LINE *Line, ELLIPSOID *Ellipsoid)
{
double vx,vy,vz,lx0,ly0,lz0,ex0,ey0,ez0,a,b,c;
register double sroot,asqr,bsqr,csqr,vxdx,vydy,vzdz;
double t;
double dx,dy,dz;
double vxsqr,vysqr,vzsqr,dxsqr,dysqr,dzsqr;
double r,s,t1,t2;
EllipsoidIntersectionAttempts++;
vx=Line->Direction.x; vy=Line->Direction.y; vz=Line->Direction.z;
lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
ex0=Ellipsoid->Centre.x; ey0=Ellipsoid->Centre.y; ez0=Ellipsoid->Centre.z;
a=Ellipsoid->Radius.x; b=Ellipsoid->Radius.y; c=Ellipsoid->Radius.z;
asqr=a*a; bsqr=b*b; csqr=c*c; /* These are used more than once, */
dx=lx0-ex0; dy=ly0-ey0; dz=lz0-ez0; /* thus much time is saved by not */
vxdx=vx*dx; vydy=vy*dy; vzdz=vz*dz; /* performing the same operation */
dxsqr=dx*dx; dysqr=dy*dy; dzsqr=dz*dz; /* twice (or even more). */
vxsqr=vx*vx; vysqr=vy*vy; vzsqr=vz*vz;
s=-((vxdx/asqr)+(vydy/bsqr)+(vzdz/csqr));
sroot=csqr*(vxsqr*bsqr+vysqr*asqr)+vzsqr*asqr*bsqr;
sroot+=2*(vxdx*(csqr*vydy+bsqr*vzdz)+asqr*vydy*vzdz);
sroot-=vxsqr*(csqr*dysqr+bsqr*dzsqr);
sroot-=vysqr*(csqr*dxsqr+asqr*dzsqr);
sroot-=vzsqr*(bsqr*dxsqr+asqr*dysqr);
if(sroot>=0) {
sroot=sqrt(sroot)/(a*b*c);
r=((vxsqr/asqr)+(vysqr/bsqr)+(vzsqr/csqr));
t1=(s+sroot)/r;
t2=(s-sroot)/r;
if((t1>EPSILON)&&((t1<t2)||(t2<EPSILON))) {
t=t1;
EllipsoidIntersections++;
}
else {
if((t2>EPSILON)&&((t2<t1)||(t1<EPSILON))) {
t=t2;
EllipsoidIntersections++;
}
else { t=(double) 0; }
}
}
else t=0.0;
return(t);
}
double Intersect_Triangle(LINE *Line, TRIANGLE *Triangle)
{
double t=0.0,ttmp,pnx,pny,pnz,a;
double lvx,lvy,lvz,lx0,ly0,lz0,angle1,angle2,angle3;
VECTOR v1,v2,v3;
POINT ip;
TriangleIntersectionAttempts++;
lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
pnx=Triangle->Plane.Normal.x; pny=Triangle->Plane.Normal.y; pnz=Triangle->Plane.Normal.z; a=Triangle->Plane.a;
ttmp=-(pnx*lx0+pny*ly0+pnz*lz0-a)/(pnx*lvx+pny*lvy+pnz*lvz);
if(ttmp>EPSILON) {
ip.x=Line->Origin.x+ttmp*Line->Direction.x; /* ip = intersection point */
ip.y=Line->Origin.y+ttmp*Line->Direction.y;
ip.z=Line->Origin.z+ttmp*Line->Direction.z;
if((ip.x>Triangle->Min.x)&&(ip.x<Triangle->Max.x)&&(ip.y>Triangle->Min.y)&&(ip.y<Triangle->Max.y)&&(ip.z>Triangle->Min.z)&&(ip.z<Triangle->Max.z)) {
CreateVector(&v1,Triangle->Corners[1].x-Triangle->Corners[0].x,Triangle->Corners[1].y-Triangle->Corners[0].y,Triangle->Corners[1].z-Triangle->Corners[0].z);
CreateVector(&v2,Triangle->Corners[2].x-Triangle->Corners[0].x,Triangle->Corners[2].y-Triangle->Corners[0].y,Triangle->Corners[2].z-Triangle->Corners[0].z);
angle1=VectorsAngle(&v1,&v2);
CreateVector(&v3,ip.x-Triangle->Corners[0].x,ip.y-Triangle->Corners[0].y,ip.z-Triangle->Corners[0].z);
angle2=VectorsAngle(&v1,&v3);
angle3=VectorsAngle(&v2,&v3);
if((angle2+angle3)<=(angle1+EPSILON)) {
CreateVector(&v1,Triangle->Corners[1].x-Triangle->Corners[2].x,Triangle->Corners[1].y-Triangle->Corners[2].y,Triangle->Corners[1].z-Triangle->Corners[2].z);
CreateVector(&v2,Triangle->Corners[0].x-Triangle->Corners[2].x,Triangle->Corners[0].y-Triangle->Corners[2].y,Triangle->Corners[0].z-Triangle->Corners[2].z);
angle1=VectorsAngle(&v1,&v2);
CreateVector(&v3,ip.x-Triangle->Corners[2].x,ip.y-Triangle->Corners[2].y,ip.z-Triangle->Corners[2].z);
angle2=VectorsAngle(&v1,&v3);
angle3=VectorsAngle(&v2,&v3);
if((angle2+angle3)<=(angle1+EPSILON)) {
t=ttmp;
TriangleIntersections++;
}
}
}
}
return(t);
}
double Intersect_Box(LINE *Line, BOX *Box)
{
double t=0.0,told=-1.0,ttmp,pnx,pny,pnz,a;
double lvx,lvy,lvz,lx0,ly0,lz0;
int i;
POINT ip;
BoxIntersectionAttempts++;
lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
for(i=0;i<6;i++) {
pnx=Box->Planes[i].Normal.x; pny=Box->Planes[i].Normal.y; pnz=Box->Planes[i].Normal.z; a=Box->Planes[i].a;
ttmp=-(pnx*lx0+pny*ly0+pnz*lz0-a)/(pnx*lvx+pny*lvy+pnz*lvz);
if((ttmp>EPSILON)&&((told<0.0)||(ttmp<told))) {
switch(i) {
case 0:
case 1:
ip.y=Line->Origin.y+ttmp*Line->Direction.y;
ip.z=Line->Origin.z+ttmp*Line->Direction.z;
if((ip.y>=Box->Corners[0].y)&&(ip.y<=Box->Corners[1].y)&&(ip.z>=Box->Corners[0].z)&&(ip.z<=Box->Corners[1].z)) {
told=ttmp;
}
break;
case 2:
case 3:
ip.x=Line->Origin.x+ttmp*Line->Direction.x;
ip.z=Line->Origin.z+ttmp*Line->Direction.z;
if((ip.x>=Box->Corners[0].x)&&(ip.x<=Box->Corners[1].x)&&(ip.z>=Box->Corners[0].z)&&(ip.z<=Box->Corners[1].z)) {
told=ttmp;
}
break;
case 4:
case 5:
ip.x=Line->Origin.x+ttmp*Line->Direction.x;
ip.y=Line->Origin.y+ttmp*Line->Direction.y;
if((ip.x>=Box->Corners[0].x)&&(ip.x<=Box->Corners[1].x)&&(ip.y>=Box->Corners[0].y)&&(ip.y<=Box->Corners[1].y)) {
told=ttmp;
}
break;
}
}
}
if(told>EPSILON) {
t=told;
BoxIntersections++;
}
return(t);
}
double Intersect_Disc(LINE *Line, DISC *Disc)
{
register double t,pnx,pny,pnz,a,tmp;
double lvx,lvy,lvz,lx0,ly0,lz0,dx,dy,dz;
DiscIntersectionAttempts++;
lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
pnx=Disc->Plane.Normal.x; pny=Disc->Plane.Normal.y; pnz=Disc->Plane.Normal.z; a=Disc->Plane.a;
tmp=pnx*lvx+pny*lvy+pnz*lvz;
if(tmp!=0.0) {
t=-(pnx*lx0+pny*ly0+pnz*lz0-a)/tmp;
if(t<EPSILON) { t=0.0; }
else {
dx=(Line->Origin.x+Line->Direction.x*t)-Disc->Centre.x;
dy=(Line->Origin.y+Line->Direction.y*t)-Disc->Centre.y;
dz=(Line->Origin.z+Line->Direction.z*t)-Disc->Centre.z;
if((dx*dx+dy*dy+dz*dz)>(Disc->r*Disc->r)) { t=0.0; }
else { DiscIntersections++; }
}
}
else { t=0.0; }
return(t);
}
double Intersect_Cylinder(LINE *Line, CYLINDER *Cylinder)
{
register double t=0.0,t1,t2,t3,t4,rsqr,lx0lvx,ly0lvy;
double lvx,lvy,lvz,lx0,ly0,lz0,tmp,tmp2,tmp3;
POINT ip1,ip2;
CylinderIntersectionAttempts++;
lvx=Line->Direction.x; lvy=Line->Direction.y; lvz=Line->Direction.z;
lx0=Line->Origin.x; ly0=Line->Origin.y; lz0=Line->Origin.z;
rsqr=Cylinder->r*Cylinder->r;
lx0lvx=lx0*lvx; ly0lvy=ly0*lvy;
tmp=2.0*lx0lvx*ly0lvy-lvx*lvx*ly0*ly0-lvy*lvy*lx0*lx0+lvx*lvx*rsqr+lvy*lvy*rsqr;
if(tmp>=0.0) {
tmp=sqrt(tmp);
tmp2=-(lx0lvx+ly0lvy);
tmp3=1.0/(lvx*lvx+lvy*lvy);
t1=(tmp2+tmp)*tmp3;
t2=(tmp2-tmp)*tmp3;
ip1.x=lx0+lvx*t1; ip1.y=ly0+lvy*t1; ip1.z=lz0+lvz*t1;
ip2.x=lx0+lvx*t2; ip2.y=ly0+lvy*t2; ip2.z=lz0+lvz*t2;
tmp=ip1.x*Cylinder->Discs[0].Plane.Normal.x+ip1.y*Cylinder->Discs[0].Plane.Normal.y+ip1.z*Cylinder->Discs[0].Plane.Normal.z;
if(tmp>Cylinder->Discs[0].Plane.a) { t1=0.0; }
else {
tmp=ip1.x*Cylinder->Discs[1].Plane.Normal.x+ip1.y*Cylinder->Discs[1].Plane.Normal.y+ip1.z*Cylinder->Discs[1].Plane.Normal.z;
if(tmp>Cylinder->Discs[1].Plane.a) { t1=0.0; }
}
tmp=ip2.x*Cylinder->Discs[0].Plane.Normal.x+ip2.y*Cylinder->Discs[0].Plane.Normal.y+ip2.z*Cylinder->Discs[0].Plane.Normal.z;
if(tmp>Cylinder->Discs[0].Plane.a) { t2=0.0; }
else {
tmp=ip2.x*Cylinder->Discs[1].Plane.Normal.x+ip2.y*Cylinder->Discs[1].Plane.Normal.y+ip2.z*Cylinder->Discs[1].Plane.Normal.z;
if(tmp>Cylinder->Discs[1].Plane.a) { t2=0.0; }
}
t=0.0;
if((t1<EPSILON)||(t2<EPSILON)) {
t3=Intersect_Disc(Line, &Cylinder->Discs[0]);
t4=Intersect_Disc(Line, &Cylinder->Discs[1]);
if(t1>EPSILON) { t=t1; }
if((t2>EPSILON)&&((t2<t)||(t<EPSILON))) { t=t2; }
if((t3>EPSILON)&&((t3<t)||(t<EPSILON))) { t=t3; }
if((t4>EPSILON)&&((t4<t)||(t<EPSILON))) { t=t4; }
}
else {
if(t1>EPSILON) { t=t1; }
if((t2>EPSILON)&&((t2<t)||(t<EPSILON))) { t=t2; }
}
if(t>EPSILON) { CylinderIntersections++; }
}
return(t);
}